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Summary 


There  has  been  increased  interest  recently  in  feed-back  methods  for  reliable,  robust, 
efficient  computational  methods  in  mechanics.  We  will  outline  the  construction  of 
such  methods  for  a  class  of  problems  describing  special  (anti-plane  shear)  deforma¬ 
tions  of  bars  of  rectangular  or  arched  cross  section.  In  particular,  we  will  show  how 
to  reduce  the  dimension  of  the  underlying  problem  “adaptively”.  For  brittle  or  linear 
materials,  this  method  is  adaptive  (optimal  in  rate  of  convergence). 

We  shall  emphasize  the  computational  aspects  that  have  practical  import  to  the 
performance  of  this  method,  such  as  the  construction  of  a  posteriori  error  estimators 
that  are  simple  to  compute,  the  selection  of  basis  functions  in  the  dimensional  reduc¬ 
tion  and  the  heuristic  principle  for  extension.  We  will  illustrate  these  concepts  with 
computations.  /  /  J  j  ( _ 


1.  Introduction 

Since  the  dawn  of  structural  mechanics,  the  scientific  and  engineering  community  has 
pursued  the  idea  of  replacing  a  full  system  of  equations  governing  the  behaviour  of 
“thin”  structures  by  ones  posed  over  their  mid  surface.  We  roughly  categorize  these 
as  follows: 

•  a  physical  modelling  approach  where  one  uses  a  priori  physical  assumptions 
along  with  appropriate  variational  principles  to  arrive  at  the  dimensionally  re¬ 
duced  equations.  See  (16]  for  a  survey. 

•  an  asymptotic  expansion  approach.  Based  on  -first-  formal  asymptotic  expan¬ 
sion  in  a  parameter  measuring  the  “thin”  dimension,  one  identifies  dimension- 
ally  reduced  equations  -  often  “justifying”  those  arrived  at  above  by  proving 
convergence  as  this  “thin”  parameter  tends  to  zero  -  under  sometimes  very 
restrictive  assumptions.  See  (4]  for  a  survey. 

•  an  energyasymptotic  approach,  which  is  the  one  we  will  pursue  here,  the  basic 
idea  is  to  find  a  minimizer  of  the  given  energy  functional  in  a  proper  (but 
still  of  infinite  dimension)  subspace  Vs  which  is  characterized  by  some  basis 
functions  {V*>}yLo: 

vN  = 
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Here  we  assume  that  the  midsurface  is  parametrized  by  the  equation  xn  —  0 
and  that  the  thickness  of  the  domain  is  d(x').  Thus  the  model  of  order  N  of 
reduced  dimension  was  introduced.  In  [6]  the  choice  of  {V>;}yt0  was  made  to 
yield  optimal  rates  of  convergence  as  d  J.  0  and  as  by-product  also  as  N  ->  oo. 
This  approach  goes  back  at  least  as  far  as  to  Poniatovskii,  (15),  and  Kantorovich, 
(12). 

Undoubtedly,  other  relevant  categorizations  may  be  thought  of.  If  one  uses  piecewise 
polynomials  for  ipj  and  discretizes  the  Galerkin  equations  posed  over  VN  fully  by  a 
Galerkin  finite  element  method,  again  using  piecewise  polynomials,  we  may  identify 
the  resulting  discretization  as  one  of  an  anisotropic  p  version  using  a  tensorial  basis 
split  along  the  “thin”  direction.  We  will  remain  free  to  choose  a  full  discretization, 
however.  For  one  dimension  there  exist  robust  nonlinear  boundary  value  problem 
solvers.  In  addition,  the  coefficient  terms  c;  do  not  coincide  with  the  terms  in  the 
formal  asymptotic  expansion  or  the  reduced  solution  in  the  first  approach  above. 

In  an  earlier  paper  (6),  the  method  of  dimensional  reduction  for  quasilinear  bound¬ 
ary  value  problems  was  introduced.  It  was  observed  that  the  coefficient  functions  Cj 
were  concentrated  at  boundary  and  interior  layers  for  large  j/d.  A  generalization  was 
proposed  which  allows  for  the  possibility  of  different  order  of  dimensionally  reduced 
models  in  different  parts  of  the  underlying  domain. 

We  shall  confine  our  study  to  the  class  of  problems  describing  special  (anti-plane 
shear)  deformations  of  bars  of  rectangular  or  arched  cross  section.  In  fig.  1  we  depict 
an  arched  cross  section  with  the  midsurface  being  described  by  the  equation  -  in 
polar  coordinates 

^mid  •  r  =  ro(^)l  0  £  (0,w) 

The  reference  configuration  for  the  infinite  bar  then  lias  a  cross  section 

=  (M)  :  r  e  ro(0)  +  (-d(0),  +d(0));  0  <=  (0,u/) 

We  briefly  derive  the  governing  equations.  During  an  anti-plane  shear  a  prism  with 
generators  parallel  to  x3  axis  undergo  the  deformation  {x\,x2,x3)  *-*  (xt,i2,is  |- 
U( xj,xj))  with  a  deformation  gradient 

1  0  0 

0  1  0 
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and  the  left  Cauchy-Green  strain  tensor:  B  .  =  / fT .  For  a  neo-Hookean,  hypere¬ 
lastic  material  -  i.e.  nonlinearly  elastic,  homogeneous,  isotropic  and  incompressible 
material  -  the  strain  energy  is  given  by  W (trace#)  =  IF(3  +  |Vt/|2)  which  in  turn 
determines  the  Piola-Kirchhoff  stress  tensor:  S  —  —pf~T  +  2Wr'(traceB)/  up  to  ar¬ 
bitrary  pressure  p  so  that  Su  =  S22  =  S33  =  2W^'(tracefl)  —  p,  S12  =  S21  =  0,  and 
S31  =  2fP'(traceJ?)|^,  S32  =  2W1  (traceS)  From  the  equilibrium  equations: 

tgdx(Su)=0 

follows 

~{2W'(3  +  |Vf/|2)~  +  ^(2^'(3  +  |V£/|2)|^  =  const. 

Identifying  F(|VI7|2)  with  2W'(3  +  |Vt/|2),  the  equations  are  then, 

divF(|Vt/|2)V£/  =  g 

with  zero  displacements  prescribed  along  0  =  0  and  0  =  u>  and  traction  (3(0)d  pre¬ 
scribed  along  r+  =  (0,u;)  x  r0  +  d  and  T+  =  (0,  uj)  x  r0  —  d.  F  is  constitutively  given 
and  we  take  simple  power  laws: 

F’(t)  =  1  +  tn,for  some  n  G  N,V<  >  0  (1.1) 


Introduce  the  scaling 


Let  u(p,0)  =  U[x).  Then  we  relist  the  equation  for  the  linear  case  F  =  1: 
r2AU  =  [(r0  +  dp)2  +  (r„  +  d'p)2Jup„ 

+{^[(r°  +  dp)  -  (r”0  +  <Fp)\  +  Jj((ro  +  d!p)d!\}up  (1.2) 
2 

+  d'p)u„,  +  u„) 

=  (ro  4-  dp)2  f 

For  n  >  1,  the  algebra  becomes  more  cumbersome  but  manageable. 

Since  we  use  a  Galerkin  approach  we  put  this  into  a  variational  form.  Let  w  =  1. 
Find  u  G  V  such  that 


Vv  G  V,  i4u(v)  =  G(t>) 

(1.3) 

Au(v)  =  J  F(|  Vjti  |!)  Vati  ■  Vjv  d  i6  ip 

(14) 

=  *('  0(9)|»(M)(1  +  i)  -  t>(-l,»)(l  -  i )]  M 
•'o 

(1.5) 

n=|o,i|x]-i,i| 

(1.6) 

To  =  ({0}  x|  —  1, 1()  U  ({l}x|  -  1, 1|) 

(1.7) 

r*  =  (o,i|  x  {i} 

(1.8) 

r.  =  (o,i|  x  {-i> 

(1.9) 

where 


4 


v  =  w^n+2(n)  =  {ve  ^^(njjBir,  =  0} 


v. 


cos  6  —  sin  0 
sin#  cos# 


,ld_  1  d 
ddp'  1  +  pddB 


) 


(1.10) 

(1.11) 


To  make  the  algebra  as  simple  as  possible,  we  worked  with  the  monomials  as 
basis  in  L2(  — 1, 1). 

This  problem  corresponds  not  only  to  finding  a  minimizer  in  V  for  the  energy  in 
anti-plane  shear  in  finite  elasticity,  [9]  and  [10]  but  also  to  the  torsion  problem  for  a 
bar,  see  [11]  and  [13];  then  U  denotes  the  Prandtl  stress  function.  See  [6].  It  is  also 
obtained  from  Maxwell’s  equations  for  the  magnetic  vector  potential  related  to  the 
magnetic  flux  density  B  when  the  magnetic  permeability  depends  on  |i?|2. 


2.  Dimensional  Reduction 

We  define  the  dimensionally  reduced  solution  of  order  N  to  be  the  solution  uN  in 
V}v  C  V  for  which 

Vv  €  VN  C  V ,  Aujv(v)  =  G(v)  (2.1) 

given  Vat  a  subspace  in  V  of  the  form 

VN  =  {v£V:  v(9,p)  =  £f=0c,(#)^(#)}  (2.2) 

The  family  of  subspaces  {Vjv}^=0  *s  characterized  by  the  choice  of  0  called  the 

basis  or  Ansatz  functions. 

In  [6]  these  basis  functions  were  selected  to  yield  the  optimal  rate  of  convergence 
of  ||u  -  uyv||//i  as  d  j  0.  We  thus  had  to  select  xf,  to  be  a  polynomial  of  degree  2j. 
For  F  in  (1.1)  depending  on  p  in  a  more  general  way,  [6]  indicated  that  the  same 
procedure  would  yield  ipj  to  be  a  nonpolynomial  solution  of  a  second  order  Sturm 
Liouville  problem.  (See  also  [Remark  3.9] [7].) 

Let  us  determine  the  basis  functions.  Setting  u  =  in  the  PDE  for  a 

circular  arch: 

—  [l  +  pd\2upp  +  -[l  +  pd\up  +  Uft  =  0 
yields  a  formal  asymptotic  expansion  defined  recursively  by 

ut+2)  +  {2#ui7X)  +  u?+1)}  +  {p2upI  +  Pup]  +  ute  }  =  0 

with  =  0  for  fc  <  0;  the  corresponding  BC  are  u£+1)  =  /36,i0  at  p  =  ±1.  Carrying 
this  proces  out,  one  obtains  that  ub)  can  be  written  as  a  polynomial  in  p  of  degree 
j  with  coefficients  that  are  functions  of  9  determined  intrinsically  in  the  proces.  The 
solution  to  the  Galerkin  equations  is  determined  by  the  span  of  the  basis  functions, 
which  then  is  the  polynomials  up  to  degree  N.  Using  the  techniques  in  [6]  it  is  possible 
to  prove  convergence  as  ef  j.  0.  Let  us  be  the  Nth  partial  sum  in  the  formal  asymptotic 
expansion  as  given  above,  e.g.  Let  D\  be  the  operator  defined  by  D\u  =  £$  mapping 
Dom(Di)  =  W2-2n+2(0, 1)  n  W0U"+2  -  L2n+2(0, 1).  We  get. 

Theorem  2.1  Let  n  6  Z+.  Let  u ,un  €  W1,00  be  bounded  there  independently  of  d. 
Let  0  €  Dom(Di).  Then  there  exists  CN  independent  of  d  such  that 

Ilu  -  ||ff»  <  C/tdiff+% 
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Regarding  convergence  of  the  formal  asymptotic  expansion  for  the  linear  case, 
we  showed  in  [7]  that  it  is  uniformly  convergent  if  d  is  sufficiently  small  and  /?  is 
sufficiently  smooth  and  compatible  -  it  had  to  be  more  than  analytic:  /3^2*^(0)  = 
0(2fc)(l)  =  0,  VA:  €  N  and  3 M  >  0  such  that  Vi  €  N  :  ||/3^||oo  <  M*  -  very  restrictive 
conditions.  For  the  partial  sum  up  to  index  N  to  exist,  one  needs  the  boundary 
conditions  above  for  k  <  N. 

Since,  for  a  given  practical  problem,  we  cannot  depend  on  d  being  sufficiently 
small  to  ensure  that  a  given  tolerance  criterion  can  be  satisfied  via  Theorems  2.1,  we 
have  considered  in  [6]  to  increase  N .  Again,  optimal  rates  in  this  scenario  (d  fixed, 
N  increasing)  were  established  in  [6].  From  the  computational  experience  in  [6]  and 
elsewhere,  it  became  clear  that  it  was  unnecessary  (read:  wasteful)  to  increase  N 
uniformly  everywhere  in  [0,1].  Rather,  there  were  clearly  defined  layers  (near  the 
boundary  and/or  rough  spots  in  the  load).  We  propose  to  increase  N  near  these 
layers  only  as  our  extension  procedure. 

Let  I  =  (0, 1)  =  U£,I„  and  I,  Hi,  =  0,  i  ^  j  Vt,j  6  jl,  m].  Let  M  =  (A,)”  j  be  an 
m-vector  of  nonnegative  integers  (TV,  =  no.  of  basis  functions  used  in  I,).  Consider 

Kv  =  {v  •'  u(0,p)  =  EyL0 Vj(0)xpj(p)  such  that  (2.3) 

N  =  ||A/!|oo,  Vj{0)  =  0  for  0  G  Uy>Wi.I,} 

a  subspace  of  VN.  Solving 

Vv  €  VM  C  VN,  Au}j(v)  =  G(v)  (2.4) 

for  u  €  Fj i/  is  the  generalized  dimensionally  reduced  Galerkin  problem. 

A  key  ingredient  in  the  selection  of  the  distribution  of  orders  M  -  the  local  a 
posteriori  estimators  which  will  identify  the  aforementioned  layers  -  will  be  developed 
in  the  following  section. 

3.  A  Posteriori  Error  Estimators 

Define  the  estimator  for  (0,1)  and  order  M  as 

=  H—ilwn)  (3.1) 

where  e  6  //^(fl)  is  the  solution  of 

/I  f  1  1  Qfi  1 

Vo  dd~tddeidei,'°a(v)~A “*(*)  *3-2) 

the  right  hand  side  being  the  residual  (Au  -  Au h){v).  Although  e  is  not  well  defined, 
and  Est(M)  are,  provided  the  following  solvability  condition  is  satisfied 

Vc  €  1  (0, 1)  :  / *  /?(0)2c(0)d  d9  =  (3.3) 

•'o 

/’ /'f( I  Vjuj,  |’)^£c 
-1  J0  Off 

However,  this  is  satisfied  (even  for  c  €  W01,Jn+2)  if 

1  €  «pafi({0;>JLo) 


(3.4) 
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cf.  (2.1)  and  (2.2).  This  condition  is  met  for  any  choice  of  basis  functions  with 
optimal  rates,  see  [6]. 

Similarly  define  the  local  error  estimator 

Esti(M)  =  [f  f  (i^r-)2  d  dOdp\  *  ,l<t<m  (3.5) 

-i  •'I,  a  dp 

As  in  [2]  we  define  upper  (lower)  error  estimator  to  mean 

||u  -  ujyIIh*  <  (>)Est 

Theorem  3.1  Let  u  and  uh  be  the  exact  and  dimensionally  reduced  solutions  (see 
(l.S)  and  (2.4))-  Then  Est  as  defined  tn  (8.1)  is  an  upper  estimator,  i.e. 

||u  -  u.v||w.  <  Est(M) 

In  the  language  of  [3],  Est  is  a  guaranteed  U-estimator  (G-estimator).  Also  Est  is 
asymptotically  exact  in  the  linear  case  (and  under  some  restrictions  in  the  nonlinear 
one),  i.e.  Est/\\u  —  uN  ||#i  tends  to  1  as  d  tends  to  zero  for  3  sufficiently  smooth.  For 
the  nonlinear  case,  d  not  small  one  gets  a  bound  like  Est(M)  <  C{d)\\u  -  uM\\  /T'¬ 
other  estimators  are  of  course  possible.  A  whole  class  of  estimators  can  be  intro¬ 
duced  from  this  residual  based  approach  as  indicated  in  [21]. 

The  localization  of  the  error  estimator  Est  as  defined  in  (3.5)  can  be  founded  on 
exponential  decay  of  the  solution  away  from  ’’vertical”  boundaries  and/or  rough  spots 
in  the  load  as  defined  in  (1.5).  The  generalized  Galerkin  problem  (2.4)  transforms 
to  a  system  of  O.  D.  E.s  whose  solution  can  be  bracketed  by  solutions  to  two  linear 
ones.  These  two  obey  a  “St.  Venant  principle”,  i.e.  decay  exponentially  away  from 
a  concentrated  load.  See  [8].  In  fact,  for  any  two  optimal  bases  -  spanning  the 
polynomials  up  to  degree  N  in  our  case  -  the  exponential  decay  rates  are  the  same. 
For  example  -  when  the  problem  is  symmetric  in  p  if  VN  =  span({l,p2}),  the 
two  first  decay  rates  (eigenvalues)  are  k.%  =  0  (corresponding  to  eigenfunction  1)  and 
k\  =  15,  the  latter  approximating  well  the  eigenvalue  of  the  full  problem  Ax  =  7r2  with 
respect  to  exponential  decay.  In  contrast,  if  one  omits  1  as  the  first  basis  function 
(nonoptimal),  Kq  =  0  =  k\.  For  N  >  2,  the  approximation  of  Ax  can  not  get  any 
worse.  For  TV  =  3,  we  already  get  /cx  =  9.9412  «  9.8696  =  ~2. 

It  is  possible  for  each  N  to  diagonalize  the  system  of  O.D.E.s  in  the  linear  case 
and  even  tabulate  in  advance.  However  this  is  not  possible  to  do  uniformly  for  all 
N  and  loses  importance  for  the  nonlinear  case.  Since  we  consider  it  important  to 
proceed  hierarchically,  we  instead  aim  at  small  bandwidth.  We  therefore  choose 
the  basisfunctions  as  integrals  of  Legendre  polynomials.  If  the  partial  sum  in  the 
asymptotic  series  Us  exists  (/?  6  Dom(D f)  and  sufficiently  small)  and  one  chooses  as 
the  basis  the  ones  from  the  asymptotic  expansion  (special  combinations  of  monomials 
up  to  degree  j  -  depending  on  /?),  one  needs  only  compute  the  very  last  coefficient 
cs  but  it  is  not  practical  to  require  this  set  of  circumstances. 

4.  Adaptivity  and  Computational  Aspects 

We  now  wish  to  define  our  feedback  extension  procedure  which  under  further  re¬ 
strictions  turn  out  to  be  adaptive,  i.e.  optimal  with  respect  to  convergence  rate  as 
N  — ♦  oo.  This  is  the  precise  sense  of  “adaptive”  due  to  BabuSka  and  Rheinboldt  [17]. 

We  first  introduce  the  heuristic  principle  which  will  guide  us  to  an  efficient  exten¬ 
sion  procedure  based  on  the  local  a  posteriori  error  estimators. 
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Heuristic  4.1  Let  the  error  associated  with  the  generalized  dimensional  reduction  be 
estimated  by 

and  the  cost  (work)  be  estimated  by 


Then  we  aim  at  achieving 


Est}(Ni)  -  Est\{Ni  -  1)  oc  1,1,) 

by  increasing  by  a  factor  7  >  1  where  the  error-cost  quotient  is  maximal. 

Reasoning:  Minimizing  the  error  at  fixed  cost  with  respect  to  Ni}  yields  via  La¬ 
grange’s  multiplier  and  a  backwards  difference  approximation  the  proportionality 
aimed  at  in  the  Heuristic. 

//// 

A  typical  choice  for  workestimate  is 

=  (ouNi  +  l)aj|I,|  (4.1) 

for  some  choice  of  positive  a,,  i=l,2. 

An  easy  computation  shows  that  adding  one  to  the  relevant  iV,  will  cost  one  order 
of  N  more  than  the  computation  of  u»  itself  and  that  if  one  multiplies  by  7,  Est  can 
be  computed  at  about  the  same  cost  as  v.y . 

Remark  4.2  Increasing  Ni  by  factor  7  >  1  where  necessary  is  sufficient  to  obtain 
a  cost  0}  computing  Est  at  the  same  order  of  magnitude  or  less  of  computing  u *  , 
however  it  is  not  sufficient  to  increase  Ni  by  adding  1. 

Introduce  the  convergence  rate  function  $  for  any  given  positive  integer  N : 

*(N)  =  |J,nf  A,  Ln/  11“  -  VlUl 

||A/||oo  =N  v€Vm 

We  can  now  state  the  following  theorem  about  adaptivity  w.r.t.  convergence  rate: 

Theorem  4.3  Let  u  and  be  the  exact  and  dimensionally  reduced  solutions  and 
have  gradients  bounded  uniformly  in  d.  Let  u  6  Then 

||u-«*||  <C*{N) 

The  anisotropic  Sobolev  space  was  introduced  in  [6]  and  the  superpair  in¬ 

dicates  up  to  what  order  generalized  partial  derivatives  (with  respect  to  the  corre¬ 
sponding  variable  in  the  subpair)  are  Lj  integrable. 

It  is  possible  to  extend  this  result  to  the  case  where  m  — *  00  as  N  — ►  00.  We 
here  envision  computing  Est{$),  V0  6  I  and  based  on  this  decide  whether  a  bisection 
and  an  increase  of  N  in  a  half  interval  is  beneficial  (this  is  possible  using  obvious 
modifications  to  Heuristic  4.1). 

From  an  implementational  point  of  view,  the  nice  mathematical  properties  of  Est 
established  in  the  previous  section  will  not  suffice,  since  finding  e  as  a  solution  of  a 
second  order  O.D.E.  might  be  too  costly.  It  is  possible  to  find  formulae  that  can  be 
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used  to  compute  e  and  Est  in  practice  using  that  e  has  to  be  polynomial  in  p  which 
is  orthogonal  to  VN  in  the  /_!j  Jq  upvfi  d  d0dp  inner  product.  Let  N  =  0,  e.g.  Then 

Est2  {6)  =  ^(32d3 

Note  that  we  selected  the  monomials  as  basis  merely  to  be  able  to  formulate  Au(u) 
easily;  a  change  of  basis  within  the  same  span  requires  a  linear  transformation  only 
in  order  to  modify  the  formulae  for  the  new  choice  of  basis  functions. 

We  next  describe  briefly  our  computational  experience  with  the  estimator  above 
for  the  case  of  a  rectangular  cross  section  and  a  nonlinearly  elastic,  brittle  material. 
Results  for  the  linear  and  nonlinear  arch  will  be  reported  elsewhere.  We  complement 
the  experiments  in  [21],  where  a  constant  load  was  treated  in  the  linear  case  and  the 
viability  of  a  similar  estimator  was  established  even  for  such  noncompatible  loads. 
We  solve  the  Galerkin  system  of  ODE’s  (recast  into  first  order  standard  form)  using 
the  program  NLTPBVP  developed  at  Univ.  of  Maryland,  College  Park  by  V.  Majer 
and  I.  BabuSka,  see  f  14]  and  [l] .  We  describe  the  performance  in  the  case  of  smooth 
compatible  (5  and  n  <  0.  The  exact  solution  being  unknown,  we  in  stead  use  a 
sufficiently  high  number  of  terms  in  the  asymptotic  series.  We  take  /?  =  A7r2(l  + 
3A27r2cos27rx)sin7ra:.  In  this  case,  the  reduced  solution  is  C0  =  Asin7ri.  Taking 
N  +  1  terms  in  the  asymptotic  expansion  for  purposes  of  computing  the  error  for  u n 
is  justifiable,  since  ||u  —  u^H^i  =  |]uasy,/v+i  —  un||h'(1  +  0[d2)).  We  used  SMP,  cf. 
[19]  to  generate  the  terms  and  FORTRAN  programs  for  these. 

First,  in  Table  1  we  see  the  results  for  varying  d  in  the  case  N  =  0:  9,  the 
effectivity  index  is  the  ratio  between  the  estimated  error  and  the  (approximated) 
true  error  measured  in  Hl.  The  estimator  is  within  20  %  for  domains  twice  as  thick 
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A 

Relative  error 
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0 

.1414 

1.273 

.8109 

1 

.1414 

.3461 

1.0232 

2 

.1414 

6.635E-02 

1.0550 

3 

.1414 

1.223E-02 

1.0583 

4 

.1414 

2.203E-03 

1.0587 

5 

.1414 

3.935E-04 

1.0587 

6 

.1414 

7.006E-05 

1.0587 

Table  1:  Eflictivity  indices  for  N  =  0. 


as  long  and  within  6%  when  the  domain  is  roughly  as  thick  as  long.  Est  is  not 
asymptotically  exact  for  this  nonlinear  problem.  This  overestimation  is  but  6%.  The 
results  for  N  =  1  yield  9  values  ranging  from  about  .44  at  d  =  1.  to  1.06  at  d  =  10~2 
and  smaller.  The  efficiency  of  Est  will  deteriorate  if  one  moves  in  any  of  the  following 
two  directions:  1)  towards  larger  values  of  /?  or  d  and/or  2)  towards  less  smooth,  less 
compatible  boundary  load  /?.  In  the  former  case,  the  examples  computed  so  far 
indicate  an  acceptable  performance  for  all  practical  values.  In  the  latter  case,  the 
localization  effect  that  we  will  describe  later  will  be  of  importance,  as  then  Est  will 
be  able  to  recover  as  one  refines  (introduces  new  subintervals  in  which  we  can  change 
the  orders  of  the  model  separately)  near  a  singularity. 
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In  conclusion,  the  energyasymptotic  method  for  dimensional  reduction  lends  itself 
well  to  adaptivity,  not  only  in  the  sense  that  the  Galerkin  system  of  ODE’s  is  solved 
using  an  adaptive  ODE  solver  but  also  in  the  sense  that  it  is  possible  and  desir¬ 
able  to  carry  varying  orders  of  models  in  differing  parts  of  the  underlying  domain. 
Reliable,  computable  a  posteriori  error  estimators  were  constructed.  The  general  ap¬ 
proach  need  not  be  conceptually  confined  to  the  class  of  problems  treated  here.  Our 
computations  suggest  a  practical  confirmation  and  viability  of  many  of  the  features 
described  here  of  this  method. 
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